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ABSTRACT 

We propose to determine the optimal softening length in A-body halo simulations by minimizing 
the ensemble-average acceleration error at a small radius ro. This strategy ensures that the error never 
exceeds the optimal value beyond ro- Furthermore, we derive semi-analytic formulae for calculating 
the acceleration error due to the discreteness of particles and softened gravity, which are validated 
by direct iV-body force calculations. We estimate that current state-of-the-art halo simulations suffer 
> 6% acceleration error at 1% of the halo virial radius. The error grows rapidly toward the center 
and could contribute significantly to the uncertainties of inner halo properties. 
Subject headings: galaxies: halos — methods: analytical — methods: TV-body simulations 



1. INTRODUCTION 

./V-body simulations use discrete particles to trace the 
phase-space evolution of a continuous density field un- 
der the influence of gravity. They have broad appli- 
cations in modern cosmology that r ange from struc- 
tures beyond galaxy clu sters (e.g. IPadilla fe Bauehl 
2002i IBahcall et all |2004) to earth-mass dark-matter 
halos emerging as the first ob jects in the universe 
Ipiemand. Moore, fc Stadelll2005|) . 

Because of their Monte Carlo nature, iV-body codes 
have to soften the gravity to subdue destructive effects 
of strong t wo-body scatte rings and to increase numerical 
efficiency (De hnerj l2001jh Although softening reduces 
the variance of t he force from discrete particles, it also 
introduces a bias l)Merrittll996y) . The bias increases with 
the softening length, while the variance the opposite. As 
such, there must exist an optimal softening length that 
reaches the best compromise between the bias error and 
variance error. 

Suitable softening lengths are often searched 
through convergence tests with A-bod y halo simu- 
lation s ( Navarro. Frenk. fc White! 119961 iMoore et alJ 
19981 lsiffinter et all Il998t iKnebe et all 1200(1 
Fnkushige fc Ma.kiuoll200lHPower et al.H2003D . Tt should 



be emphasized that a proper softening length must be 
matched with other simulation parameters such as the 
time step. For example, with a poor combination of the 
softening length and tim e step, the inner halo profile 
could become core-like ijFukushige. Kawai. fc Makinol 
12004 hereafter FKM04). It is not practical to search 
for every possible softening length and its matching 
simulation parameters using A-body simulations. Thus, 
the resulting softening length may not be optimal, and 
it is not clear what physical error bounds this softening 
le ngth imp oses. 

iMerrittJ l)1996|) devised a more efficient and objective 
method, which requires the optimal softening length to 
minimize the mean integrated square error of the force: 



MISE = j p(r)(|F(r) - F Truc (r)| 2 )dr, 



(1) 



where (...) stands for an ensemble average, and p(r) is 
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the continuous density. The true force refers to the New- 
tonian result in the continuous density field. The MISE 
is effectively a sum of mass-weighted square bias and 
variance, i.e. 



MISE= / p(r)|(F(r)) 



jiTru 



! (r)| 2 dr- 



p(r)[(F 2 (r))-(F(r)) 2 ]dr. 



For halo simulations, the bias error is significant only 
at the center of the halo, while the variance error de- 
creases relatively slowly outward. Hence, by minimizing 
the MISE one tends to allow large bias errors in the cen- 
ter in order to match small variance errors integrated 
over the whole halo, which may not be desirable for halo 
simulations. To ensure the accuracy of the density profile 
beyond a small radius ro, one may require that the bias 
and variance errors, rather than the integrated ones, be 
both less than a threshold for r > tq. This is the basis 
of our approach for optimizing the softening length. 

The MISE method has been implemented for halos 
with A-body force evaluations at all radii (e .g. iMerrittl 
ll996UAthanassoula et aTll2000tlDeTmenll2001() . which de- 
mands much less run time than full A-body halo simu- 
lations. For a very large number of particles, however, a 
direct TV-body summation of forces could still be time- 
consuming. To improve the efficiency, we derive semi- 
analytic expressions for the ensemble-average bias and 
variance errors assuming a Poisson sampling of the halo 
density profile with particles. In this work, we use spher- 
ically symmetric halos as targets for calculating the ac- 
celeration error and optimizing the softening length. Our 
method may be generalized for broader applications. 

2. ACCELERATION BIAS 

For convenience, we consider the acceleration bias and 
variance and express lengths in units of the halo virial 
radius r v . We set the mass of the halo within r v to 1 , so 
that the mass of each particle is m p — A -1 , where A is 
the number of particles within r v . 

The gravitational attraction between two particles can 
be generalized as F — m^f^r, e), where we have dropped 
Newton's constant. For Newtonian gravity /(r, e) = 
r~ 2 , and for Plummer softening f(r,e) = (r 2 + e 2 ) -1 . 
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Fig. 1. — The configuration for calculating the acceleration on a 
particle at a distance r from the center of the halo. 




We a lso utilize the S2 softening ijHocknev fe Eastwoodl 
1981), which treats particles as spheres of radius e/2 
with density decreasing linearly from the maximum at 
the center to zero at e/2. The S2 softening is often used 
in particle-particle parti cle-mesh codes (e.g. iCouchmanl 
199l| Uing fc FanJl994j) . 

Suppose that a particle is located at a distance r from 
the center of the halo as illustrated in Fig. ^ The mean 
radial acceleration of the particle is 



(°r) = ( / f(d, e)cos6>dm), 



p(R) ~ m p (n(R))|R|=R, 



(2) 



where dm = m p n(R)dV with n(R) and dV being the 
particle number density and volume element at R, re- 
spectively. By definition, we have 

(3) 

where p(R) is the spherically symmetric halo density. 
The number of particles n(Il)dV within the volume dV 
follows approximately the Poisson distribution so that 
its mean and variance are both N p(R)dV . In addition, 
n{J\)dV at two locations are uncorrelated, so that 

ml(n(R)n(TL'))dVdV'~p(R)p(R')dVdV' + (4) 

(5 D (R - H')N- 1 p(R)dVdV, 

where ^ D (R) is the Dirac Delta function. 
Combining equations © and J2J one gets 



7T rdf) 

sin26>d6> / f(d,e)p(R)d 2 dd, 

JO 

- - s -2dr cos 9, and the integral limit d, 
2der cos 9 + r 2 = R 2 1S , X . For Newton 
r, while for softened gravity i 



(5) 



-r 



(Or) = 71 

where R 2 = d 2 - 
solved from d 2 6 
gravity i? max = 
plus the range that softening is in effect. The bias is then 

6a = (a r ) - aj mc . (6) 

Fig. [5] shows a few examples of (a r ) under soft- 
ened gravity. The accelerations are evaluated with 
direct iV-body force summations over 10,000 ran- 
dom realizations of a Nava rro-Frenk- White (NFW, 
iNavarro. Frenk. fc Whitelll997l hereafter NFW97) halo, 
which has a concentration number c = 5. As expected, 
the acceleration bias is significant only at r less than a 
few e and is small at larger radii even where particles can 



Fig. 2. — Radial accelerations in an NFW halo of concentration 
c = 5. The halo is populated with identical particles of mass 10 -6 
out to twice its virial radius r v , and the total mass within r v is 
1. The results are obtained over 10,000 random realizations of the 
halo in each case. The true value aj*™ (solid line) is calculated for 
a continuous halo density with Newtonian gravity. Solid symbols 
correspond to the S2 softening with softening length e = 0.001 (cir- 
cles), e = 0.005 (squares), and e = 0.01 (triangles). Open symbols 
correspond to Plummer softening with e = 0.0002 (circles) and 
e = 0.001 (squares). The a T ocr behavior of the radial acceleration 
in a constant density core is drawn in a dashed line for comparison. 
The downward arrows mark the radii where the mean interparticle 
distance is equal to the softening length in the S2 cases. 



be, on average, closer than the softening length in the S2 
case. 

With softening len gths of a few percent of the virial ra- 
dius 1 (e.g. NFW97: INavarro et alJ 12004 hereafter N04, 
see Tabled, the radial acceleration at r < 0.01r v be- 
haves as if there was a constant density core, even though 
the underlying density has a central cusp with a logarith- 
mic slope of —1. Plummer softening performs worse than 
the S2 softening at the same quoted softening length, be- 
cause it has a broader softening kernel than S2 does. 

For better comparisons, we present fractional bias er- 
rors l&al/a^'" in Fig. which includes N — 10 7 re- 
sults from 10,000 realizations of the same halo. There 
is a good agreement between direct ./V-body force sum- 
mations (open symbols) and the semi-analytic results of 
equation © (solid lines). Note that the bias error does 
not depend on the number of particles within the virial 
radius. 

3. ACCELERATION VARIANCE 

At a fixed position, the acceleration fluctuates from one 
realization to another because of the discrete sampling 
of the halo by particles. These fluctuations imposes a 
sample variance error on particles' acceleration, and they 
depend on both the number of particles and softening. 
With the help of equation J3J we find the acceleration 
variance 



a" - a 



2 = ^ / siu««W 



f(d,e)p(R)d 2 dd, 
(7) 



where dg + 2dercos9 + r 2 — R 2 ut , and i? cu t formally 
extends to infinity. Since the result within the virial ra- 

1 The softening splines in NFW97 and N04 differ from S2, but 
we do not expect such differences to alter the result qualitatively. 
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Fig. 3. — Fractional acceleration bias error ( 1 6 a | /aJF 1110 , open symbols) and variance error (a^/aj^^ , solid symbols) as functions of the 
distance from the center of the halo. Solid lines that trace the symbols are the corresponding semi-analytic results using equation JHJ for 
the bias error and equation J7J for the variance error. When the interparticle distance is much larger than the softening length, one cannot 
sample the long tail error distribution well enough, even with 10,000 realizations, to measure the variance accurately. This leads to small 
discrepancies in the variance for small e at r ~ 1, which can be reduced with more realizations. The bias error does not depend on the 
number of particles N within the virial radius, and the variance error scales as TV -1 / 2 for the same softening. 



dius converges very quickly for R cut > a few r v , we set 
Rcut = 2r v to be consistent with TV-body calculations. 

We calculate fractional variance errors a a /aj ruc for the 
same set of configurations as for fractional bias errors in 
The variance errors of direct iV-body force summa- 
tions are plotted in Fig. with solid symbols. One sees 
that the TV-body results are very well traced by equa- 
tion J7J in solid lines. For a fixed N, a smaller soften- 
ing length results in larger variance errors, while, for a 
fixed softening length, the variance error scales as TV -1 / 2 . 
Plummer softening has smaller variance errors than the 
S2 softening with the same e and N because of its broader 
softening kernel. 

To access the effect of the halo profile, we calculate 
acceleration er rors using different concentration num- 
bers and using [Moore et aTl l)1999t hereafter M99) pro- 
file, which has a stronger cusp of logarithmic slope —1.5. 
The bias error does not change much with the profile. 
On the other hand, the variance error is reduced by a 
factor of 2 to 3 by boosting the concentration from c = 5 
to c = 20. The variance error of M99 halos is a factor of 
1.6 smaller than (roughly the same as) its corresponding 
NFW halos at r = 0.01 (r ~ 1), when the concentra- 
tion number is adjusted to c(M99) = [c(NFW)/1.7] 9 
(peacock fc SmitrJl2000|) . 

4. OPTIMAL SOFTENING 

As seen in Fig. [31 the bias error increases much faster 
toward the center of t he halo than t he variance error. 
Minimizing the MISE ( M errittJ |1996fl may not be op- 
timal for studying inner halo properties, because the 
MISE method prefers a relatively large softening length 
so that small variance errors spread over the entire halo 
are matched by large bias errors confined in the center. 

We propose to optimize the softening length by mini- 
mizing the mean acceleration error £ a = \Jb\ + cr 2 at a 
small radius r . In this way, one ensures that £ a never 
exceeds the optimal value at larger radii. Conversely, 
one could set an error budget at ro and determine the 
number of particles needed for a particular softening. 



TABLE 1 

Error Estimates for Selected Simulations 
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1% of the virial radius. For spline softening (NFW97, 
M99, and N04), e equals the scale beyond which the 
gravity is Newtonian. Whereas for Plummer soften- 
ing (FKM04), e is their reported softening length. 
The number of particles N within r v is approximate. 



Fig. 0| illustrates the mean acceleration error at 1% of 
the halo virial radius (dotted lines) as a function of the 
softening length for different number of particles. For an 
NFW (c = 5) halo with S2 softening, we obtain, from 
the minimum of each error curve, e op t = 0.1 liV -020 and 
£ a = 14N~ - 38 . For the same halo but with Plummer 
softening e opt = 0.0657V- 26 and £ a = 13N-° m . The 
optimal softening lengths and minimum acceleration er- 
rors of corresponding M99 halos have nearly the same 
./V- dependence as those of NFW halos but with 15-30% 
smaller prefactors, because they have smaller variance 
errors at the same soften length. 

The TV-dependence of our results are cons i stent with 
those of MISE results in lAthanassoula et al.l l|2000|) and 
iDehnenl (J2001) , despite that our halos differ from theirs 
and that we prefer smaller softening lengths. For a closer 
c omparison, we ca lculate the optimal softening length for 
a lHernauistJ l(l99 0) sphere of 10 5 particles with Plummer 
softening. We find e nnt . = 0.0049 and £ a = 30% with 
r a = 0.01. Whereas lDehnenl (|2001^ obtained e opt = 0.016 
and a mass-weighted average error of 6.7%, which lead 
to £ a = 55% at r = 0.01. With our e opt the average 
error is 7.2%. This shows that our strategy is to trade 
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slightly larger (yet tolerable) errors at large radii with 
smaller errors near the center (though 10 5 particles do 
not seem sufficient for studying inner halo properties no 
matter what softening length is used). 

5. DISCUSSION AND CONCLUSIONS 

The optimal softening lengths in this work are not truly 
optimal in the sense that they slightly depend on the as- 
sumed halo profile and that, in turn, the simulated halo 
profile could be affected by the softening length. Never- 
theless, our strategy provides guidance for choosing the 
softening length and for evaluating acceleration errors of 
simulated halos. For example, we list in Table ^ error 
estimates for four sets of halo simulations ranging from 
one of the first claims of the universal profile (NFW97, 
N < 10 4 ) to the latest investigation (FKM04, N > 10 7 ). 
The errors are typically > 6% at r = 0.01r v . Since 
the mass within 0.01x v is a few thousandths of the halo 
virial mass, even a few percent acceleration errors may 
contribute significantly to the uncertainties in the inner 
slope of the halo. In fact, N04 achieve convergent results 
only at r > 0.01r v . Hence, it may not be reliable to 
e xtrapolate to ever sm aller radii and infer a central cusp. 

IPower et all l|2003j) propose that the optimal softening 
length should satisfy the condition that the maximum 
stochastic acceleration (~ 1/Ne^ pt ) be several factors 
smaller than the mean field acceleration at r v (~ 
Roughly speaking, the acceleration variance arises from 
two sources: the stochastic acceleration as defined by 
IPower et alJ l|2003|) and Poisson fluctuations of particles 
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in the halo. Since the latter increases toward the cen- 
ter, setting a small stochastic acceleration at r v does not 
always guarantee small errors near the center. 

Our criterion for the optimal softening length imposes 
a strict upper limit on the mean acceleration error at 
r > tq. Direct error bounds on halo properties may be 
more useful for interpreting simulation results, but, to 
identify the most accurate simulation upon which error 
estimates of other simulations can be based, one needs a 
gauge like the mean acceleration error. 

To optimize the softening length for general density 
fields or halos with sub-structures or asymmetries, one 
must derive equations (0 and (0 without assuming 
spherical symmetry for the density. Moreover, one should 
generally minimize the acceleration error in high density 
regions where the error tends to be large. 

From NFW97 to N04, the acceleration error at 0.01r v 
is reduced from 42% to 8.6%, and one starts to see a 
continuously changing (logarithmic) central density slope 
instead of a constant slope of —1 in an NFW halo. Hence, 
it will be interesting to see how the results evolve when 
the error is further reduced to well below 1% with N ^> 
10 s . A recent investigation on two-body relaxation also 
suggests that 3> 10 8 particles are req uired to fait hfully 
model the very inner regions of halos l|El-Zantll2005|) . 
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